R Notebook

library(readr)
library(readxl)
library(tidyverse)
library(ggplot2)
library(hrbrthemes)
library(dplyr)
library(tidyr)
library(viridis)
library(ggpubr)
setwd("~/PostDoc/human-metnet/results/results_analysis_in_R_rmd")
dflist_names = c('fold_change_aritmetic', 'fold_change_geometric', 'fold_change_quadratic', 'fold_change_harmonic',
    'delta_aritmetic', 'delta_geometric', 'delta_quadratic', 'delta_harmonic')

read_sheets <- function(a_sheet){df <- read_excel('delta_fc_centralidades.xlsx', sheet = a_sheet) %>% column_to_rownames("...1")
                  return(df)}

map(dflist_names, read_sheets) %>% purrr::set_names(dflist_names) -> list_of_cents
names(list_of_cents)

list_of_cents$fold_change_aritmetic %>% names

Density plot of centralities

plot_density <- function(df, titulo ){cut = 0.01
                             df %>% gather() -> df_gathered
                              df_gathered['value'] %>% filter(abs(value) > cut) %>%
                              ggplot( aes(x=value)) +
                                    geom_density(fill="deepskyblue4", color="azure4", alpha=.8)+
                                       theme(legend.position="top") +
                                    ylab("Density") + ggtitle(titulo) +
                                    xlab("Centrality") + geom_vline(xintercept=0, linetype="dashed", color = "red",  size=.4) -> Density_plot_of_centralities
                              return(Density_plot_of_centralities)}



map2(list_of_cents, names(list_of_cents), plot_density)

$fold_change_aritmetic

$fold_change_geometric

$fold_change_quadratic

$fold_change_harmonic

$delta_aritmetic

$delta_geometric

$delta_quadratic

$delta_harmonic

#calcular std, kurtonis y skewness.

#centralities <-  Results[,!is.na(as.numeric(colnames(Results)))] %>% gather()

#Density_plot_of_centralities.panel <- ggarrange( NULL,Density_plot_of_centralities, nrow = 2,  heights= c(1, .4), labels = c('a','b'),
#                                                 hjust = -0.1, vjust = .8)

#Density_plot_of_centralities.panel

fba <-  read.csv('/home/alejandro/NAMU_in_progress/Figures_creation/Code/Results_Acevedo_et_al_2021/01_phpps_bipartite_fluxes_sensis_optimality/FBA_results.csv') %>% 
        select(c('ID','Name', 'Reaction', 'Flux', 'Sensitivity')) %>% column_to_rownames('ID')



append_fba <- function(cent){
                                hola <- rownames_to_column(fba)  %>% inner_join(rownames_to_column(cent)) %>% column_to_rownames('rowname')
                                return(hola)
}



map(list_of_cents, append_fba) -> fba_list_of_cents

Joining, by = "rowname"
Joining, by = "rowname"
Joining, by = "rowname"
Joining, by = "rowname"
Joining, by = "rowname"
Joining, by = "rowname"
Joining, by = "rowname"
Joining, by = "rowname"

Hierarchical clustering

library(tidyverse)
library(scales)


assign_node_type <- function(df){
data <- df  %>% rownames_to_column('ID') %>% filter(!str_detect(ID, regex('AA_|AA2', ignore_case = F)))  %>%
         mutate(`Node type` = ifelse(str_detect(Reaction, regex('\\[[a-z]A\\]',              ignore_case = T)), 'Astrocyte', NA)) %>% 
         mutate(`Node type` = ifelse(str_detect(Reaction, regex('\\[[a-z]N\\]',              ignore_case = T)), 'Neuron', `Node type`))    %>% 
         mutate(`Node type` = ifelse( str_detect(Reaction, regex('\\[[a-z]A\\]|\\[[a-z]N\\]',ignore_case = T), negate = T), 'Exchange', `Node type`)) %>%
         mutate(`Node type` = ifelse(str_detect(Reaction, regex('\\[[a-z]A\\]',ignore_case = T)) & 
                                       str_detect(Name, regex('DM_|Demand|sink',ignore_case = T)), 'Sink/Demand Astro', `Node type`))%>%
         mutate(`Node type` = ifelse(str_detect(Reaction, regex('\\[[a-z]N\\]',ignore_case = T)) & 
                                       str_detect(Name, regex('DM_|Demand|sink',ignore_case = T)), 'Sink/Demand Neuron', `Node type`))
                                    return(data)}



map(fba_list_of_cents, assign_node_type) -> data


clean_centralities <- function(df){ df %>% select(-c ("ID", "Name", "Reaction",  "Flux", "Sensitivity","Node type")  )}


map(data, clean_centralities) -> only_centralities

glimpse(data)

data

$fold_change_aritmetic

$fold_change_geometric

$fold_change_quadratic

$fold_change_harmonic

$delta_aritmetic

$delta_geometric

$delta_quadratic

$delta_harmonic

NA

 

get_total_centralities <- function(df){df %>% abs %>% rowSums }


map(only_centralities, get_total_centralities) -> total.abs.centrality 

get_summaryzed.optimality <- function(df){ pseudoLog10 <- function(x) { asinh(x/2)/log(10) }
                        
                                      df$Flux   %>% abs %>% as.matrix()  %>% pseudoLog10  %>% rescale(c(0,1))  -> fluxes
                                    df$Sensitivity  %>% abs %>% as.matrix()   %>% pseudoLog10  %>% rescale(c(0,1))-> sensitivities 
                                    
                                    fluxes + sensitivities  -> summaryzed.optimality
                                    return(summaryzed.optimality)
                                      
}


map( data,get_summaryzed.optimality) -> summaryzed.optimality

heatmap

library(ComplexHeatmap)

library(WGCNA)


get_correlation_matrices <-function(df){
                                 #scale_rows <- function(x){t(scale(t(x)))}
                                 df %>% transposeBigData   %>% cor(method = "kendall") -> hola
                                 return(hola)}

map(only_centralities, get_correlation_matrices)  -> my.corr.matrices



#all.equal(data$fold_change_aritmetic$`Node type` , data$delta_quadratic$`Node type`)
#all.equal(data$fold_change_harmonic$`Node type` , data$delta_geometric$`Node type`)

node_types <- data$fold_change_aritmetic$`Node type`

set.seed(1)

my.corr.matrices$fold_change_geometric -> my.corr.mat

left_annotation <- rowAnnotation(`Nodes` = node_types, col = list( `Nodes` = c(
                                    "Astrocyte"        = "deepskyblue3", 
                                    'Exchange'   = 'darkgreen',
                                    'Sink/Demand Astro'='darkgoldenrod',
                                    "Neuron"           = "brown2",
                                    'Sink/Demand Neuron'= 'blueviolet' )))

total.abs.centrality$fold_change_aritmetic  -> a_total.abs.cent
summaryzed.optimality$fold_change_aritmetic -> a_summaryzed.opti

right_annotation <-   
  rowAnnotation(gap = unit(12, "points"),Ce      = anno_barplot(bar_width = 0.01,width = unit(1.5, "cm"), border = T,a_total.abs.cent, gp = gpar(col = 'azure4')),
            Op        = anno_barplot(bar_width = 0.01,width = unit(1.5, "cm") ,border = T, a_summaryzed.opti, gp = gpar(col = 'azure4')))


ht <- Heatmap(my.corr.mat, name = "Correlation",  left_annotation =left_annotation,
                                            right_annotation=right_annotation,
                                 clustering_distance_columns  = function(m)   dist(m, method = 'euclidean'),
                                 cluster_columns              = function(x) fastcluster::hclust(dist(x), "median"),
              
                                clustering_distance_rows   = function(m)   dist(m, method = 'euclidean'),
                                 cluster_rows               = function(x) fastcluster::hclust(dist(x), "median"),
                                row_km = 2,
                                column_km = 2,
                                # row_split = paste0("cluster ", pa$clustering),
                                # column_split = paste0("cluster ", pa$clustering),
                                 border = TRUE,
                                 row_dend_width    = unit(3, "cm"),
                                 row_gap = unit(2, "mm"),
                                 column_gap = unit(2, "mm"),
                                  width = unit(10, "cm"), 
                                height = unit(10, "cm"),
                                column_title = c("Astrocytic cluster", "Neuronal cluster"),
                                column_title_gp = gpar(fontsize = 10),
                                  row_title_rot  = 0,
                                 show_column_names    = F,
                                show_row_names    = F,
                                 row_names_gp = gpar(fontsize = 8),
                                  row_title = c("Astrocytic\n cluster", "Neuronal\n cluster"),
                             
                                 row_title_gp = gpar(fontsize = 10))
my_heatmap = grid.grabExpr(draw(ht))

ht

 upper_panel <- ggarrange( Density_plot_of_centralities.panel,  my_heatmap, ncol = 2, widths = c(.4, 1), heights = c(1,.3), labels = c('','c'),
                                                 hjust = -4, vjust = .8)

upper_panel

Distribution of pairwise correlations

ggplot(quads, aes(x=Comparison, y=`Node correlation`, fill=Comparison)) + 
    geom_violin(trim=F, colour = "azure4")+  
    stat_compare_means( vjust= -1., hjust= 0, comparisons = my_comparisons, method = "wilcox.test", p.adjust.method = "bonferroni", label = "p.signif")+
   scale_fill_manual(values = alpha(c("deepskyblue3", "blueviolet", "brown2"), .8)) + 
   theme(plot.title = element_text(hjust = 0.5),legend.position="none",  axis.title.x=element_blank(),
         axis.text.y = element_text(angle = 45, hjust = 1), axis.text.x = element_text(angle = 50, hjust = 1))  + 
  #ggtitle("Neuronal self-correlations") + 
  scale_y_continuous(limits=c(-0.7, 2)) + geom_hline(yintercept=0, linetype="dashed", color = "darkred",  size=1)-> corr_comparisons

corr_comparisons

library(ProjectionBasedClustering)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(PCAtools)

Loading required package: ggrepel

Attaching package: ‘PCAtools’

The following objects are masked from ‘package:stats’:

    biplot, screeplot

library(magrittr)

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract

p <- pca(my.corr.mat)

cbind(p$rotated$PC1, p$rotated$PC2) %>% as.matrix %>% as.data.frame-> my_pca




my_pca %>%  cbind( data$`Node type`) %>% set_colnames(c('PC1','PC2','node')) -> to.pca.scatter

summaryzed.optimality -> Op
total.abs.centrality -> Centrality
Centrality -> Ce

b <- ggplot(to.pca.scatter, aes(x = PC1, y = PC2)) 
b +   scale_color_manual(labels = c("Astro", "Exch",'Neu','Si/De Ast', 'Si/De Neu'), values =  c("deepskyblue3",'darkgreen','brown2','darkgoldenrod','blueviolet')) + 
  theme(legend.position="right", legend.box = "vertical")+ geom_point(aes(size = Op, color = node))-> pca_node_RedCosts




b <- ggplot(to.pca.scatter, aes(x = PC1, y = PC2)) 
b +   scale_color_manual(labels = c("Astro", "Exch",'Neu','Si/De Ast', 'Si/De Neu'), values =  c("deepskyblue3",'darkgreen','brown2','darkgoldenrod','blueviolet')) +
  theme(legend.position="right", legend.box = "vertical")+ geom_point(aes( size =  Ce,color = node))-> pca_node_centrality

bottom_panel <- ggarrange( corr_comparisons,pca_node_centrality, pca_node_RedCosts  , ncol = 3, widths  = c(.4,1,1) , labels = c('d','e','f'),
                                                 hjust = 0, vjust = 1)

Removed 3 rows containing non-finite values (stat_ydensity).Removed 3 rows containing non-finite values (stat_signif).Removed 26 rows containing missing values (geom_violin).

ggsave(file="panel_without_graph.png", plot=panel_without_graph, width=13, height=9, dpi = 320)

Node types for Networkx (python) attributes

library(tidyverse)
library(magrittr)

Results <- read_csv("Results.csv")

data_Networkx <- Results%>% 
         mutate(`Node_type` = ifelse(str_detect(Formula, regex('\\[[a-z]A\\]',              ignore_case = T)), 'Astrocyte', NA)) %>% 
         mutate(`Node_type` = ifelse(str_detect(Formula, regex('\\[[a-z]N\\]',              ignore_case = T)), 'Neuron', `Node_type`))    %>% 
         mutate(`Node_type` = ifelse( str_detect(Formula, regex('\\[[a-z]A\\]|\\[[a-z]N\\]',ignore_case = T), negate = T), 'Exchange', `Node_type`)) %>%
         mutate(`Node_type` = ifelse(str_detect(Formula, regex('\\[[a-z]A\\]',ignore_case = T)) & 
                                       str_detect(Name, regex('DM_|Demand|sink',ignore_case = T)), 'Sink/Demand Astro', `Node_type`))%>%
         mutate(`Node_type` = ifelse(str_detect(Formula, regex('\\[[a-z]N\\]',ignore_case = T)) & 
                                       str_detect(Name, regex('DM_|Demand|sink',ignore_case = T)), 'Sink/Demand Neuron', `Node_type`))


data_Networkx %<>% select(c(ID, Node_type))

write_csv(data_Networkx, 'data_Networkx.csv')